simCohort = function(parameters, lambda, 
										 populationData, 
                      deathFile,
                      DiffCancer
	 )	
 {

# parameters is a list with elements
# geno env geno:env size probGeno probEnv genoMiss envMiss oneEnvPerCommunity
# cancerMissRate falseCancerRate sdGroup, sdPerson, Enrollmentprobs, 


if(!all(c("geno", "env", "geno:env", "size", "probGeno", "probEnv", "genoMiss", 
			"envMiss", "oneEnvPerCommunity", "cancerMissRate", "falseCancerRate", 
			"sdGroup", "sdPerson", "Enrollmentprobs","Ncommunity") %in% 
		names(parameters)))
warning("some parameters missing")

#coefs = c(parameters[["geno"]], parameters[["env"]], parameters[["geno:env"]])

thecommunities = rownames(populationData$all)[sort(sample(dim(populationData$F)[1],
parameters$Ncommunity))]


# subset the communities
for(D in 1:length(populationData))
	populationData[[D]] = populationData[[D]][thecommunities,]
      
# find number of subjects per CSD
if (parameters$EqualCommunity==T) 
  thesample=rep(thecommunities, parameters$size/parameters$Ncommunity) 

if (parameters$EqualCommunity==F)  {

#{community = rownames(populationData$all)
population=populationData$all[,"total"]
thesample= sample(thecommunities, parameters$size, replace=T, prob=population/sum(population))
}

thepeople=data.frame(total=as.matrix(table(thesample)))

thepeople[["F"]] = rbinom(dim(thepeople)[1], thepeople$total, 
		populationData$all[rownames(thepeople), "propFemale"])
		
	thepeople$M = thepeople$total - thepeople$F 
                 
datamat = data.frame(CensusDivision = rep(rownames(thepeople), thepeople$total),
	Gender = rep( rep(c("F", "M"), length(thecommunities)), 
					c(t(thepeople[,c("F", "M")])))) 

#simulate age
theages = as.integer(colnames(populationData$F))

# assign an age group



ages = NULL
for(Dcommunity in thecommunities) {
	for(Dsex in c("F", "M")) {
		 	ages =  c(ages,  sample(theages, 
			 		thepeople[Dcommunity, Dsex],
		 		prob=populationData[[Dsex]][Dcommunity,],replace=T) )
		}
}




# ages are uniformly distributed within an age group
datamat$Age = ages + runif(parameters$size, 0, 5)

#simulate enrollment date
Enrollmentcdf=cumsum(parameters$Enrollment)/sum(parameters$Enrollment)
datamat$Enrollment=approx(c(0,Enrollmentcdf), 0:4, runif(parameters$size))$y

# assign gene and environment groups
datamat$geno = rbinom(parameters$size, 1,parameters$probGeno)
if(parameters$oneEnvPerCommunity) {
	env = rbinom(length(thecommunities), 1, parameters$probEnv)
	datamat$env=env[datamat$CensusDivision]
} else {
	datamat$env = rbinom(parameters$size, 1,parameters$probEnv) 	
}


# individual and group random effects
personRandomEffect = rnorm(parameters$size, 0, parameters$sdPerson)
groupRandomEffect = rnorm(length(thecommunities), 0, parameters$sdGroup)

datamat$rr = personRandomEffect + groupRandomEffect[datamat$CensusDivision] + 
      parameters$geno * datamat$geno +
      parameters$env * datamat$env + 
      parameters[["geno:env"]] * datamat$geno *datamat$env -
      # subtract off one half the variance to correct the mean of the log-normal
      0.5*(parameters$sdPerson^2 + parameters$sdGroup^2)
datamat$rr = exp(datamat$rr)


thefemales = datamat$Gender== "F"

dataFemales = datamat[thefemales,]
dataMales = datamat[!thefemales,]

# simulate non-cancer deaths

dataFemales$nonCancerDeath = 
  getDeath(deathFile$F, 
    dataFemales$Age, dataSim=dataFemales)

dataMales$nonCancerDeath =
  getDeath(deathFile$M, 
    dataMales$Age, dataSim=dataMales)
    


dataMales$CancerIncidence = simCancer(dataMales$rr, dataMales$Age,
   dataMales$nonCancerDeath, lambda$M)
   

dataFemales$CancerIncidence = simCancer(dataFemales$rr, dataFemales$Age,
   dataFemales$nonCancerDeath, lambda$F)

datamat = rbind(dataMales, dataFemales)



noCancer = datamat$CancerIncidence>199

# misclassify cancers	
missedCancer = rbinom(parameters$size, 1, parameters$cancerMissRate) & !noCancer
falseCancer = rbinom(parameters$size, 1, parameters$falseCancerRate) & noCancer

datamat$CancerIncidence[missedCancer] = 200
# missed cancers have uniform age distributions.
datamat$CancerIncidence[falseCancer] = runif(sum(falseCancer),  
	datamat$Age[falseCancer], datamat$nonCancerDeath[falseCancer])


	
# missclassify env and geno
whichGenoMiss = as.logical(rbinom(parameters$size, 1, parameters$genoMiss))
datamat$geno[whichGenoMiss] = !datamat$geno[whichGenoMiss] 

if(parameters$oneEnvPerCommunity) {
	whichEnvMiss = rbinom(length(thecommunities), 1, parameters$envMiss)
	env[whichEnvMiss] = !env[whichEnvMiss]
	datamat$env=env[datamat$CensusDivision]
} else {
	whichEnvMiss = as.logical(rbinom(parameters$size, 1, parameters$envMiss) )
	datamat$env[whichEnvMiss] = !datamat$env[whichEnvMiss] 
}

if(!is.null(DiffCancer)){

maxAge = max(as.integer(rownames(DiffCancer[[1]])))

# simulate cancer site
datamat$CancerIncidence[datamat$CancerIncidence > 199] = NA
datamat$CancerIncidenceAgeRound = 5*floor(pmin(maxAge, datamat$CancerIncidence)/5)


theorder = order(datamat$Gender, datamat$CancerIncidenceAgeRound)
datamat = datamat[theorder,]

cnumbers = list()
for(Dsex in c("F", "M")) {
	 cnumbers[[Dsex]] = table(datamat[datamat$Gender==Dsex, "CancerIncidenceAgeRound"])
}	 


	

CancerType = NULL

for(Dsex in c("F", "M")) {
	CancerNames = colnames(DiffCancer[[Dsex]]) [ !colnames(DiffCancer[[Dsex]])== "Total" ]
	DDiffCancer=DiffCancer[[Dsex]][, !colnames(DiffCancer[[Dsex]])== "Total" ]
	for(Dage in names(cnumbers[[Dsex]])) {
		 CancerType =  c(CancerType, 
		 		sample(CancerNames, cnumbers[[Dsex]][Dage], replace=T,
			 		prob = DDiffCancer[Dage,]) 
			)	  
	}
	CancerType = c(CancerType, 
			rep(NA, sum(is.na(datamat[datamat$Gender==Dsex, "CancerIncidenceAgeRound"])))
		)
}
datamat$CancerType = CancerType

} else {
 datamat$CancerType=1
}


# compute observed events
datamat$CancerIncidence[is.na(datamat$CancerIncidence)] = 200
datamat$Event=datamat$nonCancerDeath
datamat$Cancer= datamat$CancerIncidence < datamat$nonCancerDeath
datamat$Event[datamat$Cancer]= datamat$CancerIncidence[datamat$Cancer]

datamat
				
									
}

simCancer = function(Ri, start, death, lambda) {
# simulate cancer incidence
# lambda is a list with x being the age sequence and y being the rates
# the other arguments are vectors, one element per individual
Nindiv =  length(Ri)
result = rep(0, Nindiv)

      
# add age zero with zero population, 
# and age 200 with zero population
#  in case an individual has 
# age below the lowest cutoff                          
ageSeq = c(0, lambda$x, 200)
lambdaSeq = c(0,lambda$y)                                
Nlambda = length(lambdaSeq)
diffLambda = diff(ageSeq)

# these will contain the cumulative distribution for each individual
thisAgeSeq = thisLambda = rep(0, Nlambda)

if(!is.loaded("simCancer")) 
  dyn.load(paste("../src/simCancer", .Platform$dynlib.ext, sep=""))
result = .C("simCancer" , as.double(Ri), as.double(start), 
  as.double(death), result=as.double(result), 
  as.integer(Nindiv), as.double(ageSeq), as.double(lambdaSeq),
  as.double(diffLambda), as.double(thisAgeSeq), as.double(thisLambda),
  as.integer(Nlambda))$result
  

  result
}

getDeath = function(mortalityFile, age, dataSim) {
# simulate non-cancer death

mydata = mortalityFile        

library(survival)

mort.out<-survreg(Surv(time2,aalive)~1,data=mydata,weights=weights,dist=c("weibull"))
shape<-mort.out$coefficients
scale<-mort.out$scale


# simulate non-cancer death
nonCancerDeath = rweibull(dim(dataSim)[1], scale=exp(shape), shape=1/scale)
smaller = nonCancerDeath < age
while(any(smaller)) {
 nonCancerDeath[smaller] = 
   rweibull(sum(smaller), scale=exp(shape), shape=1/scale)
 smaller = nonCancerDeath < age
}

return(nonCancerDeath )
}


####fuction of calculating different cancer lambda

lambdaDiffCancer=function(myfile,agegroup)   {
require(foreign)
thefile = read.table(file=myfile, sep=",", header=T, as.is=T, quote="\"")

# get rid of rows with totals in them
thefile = thefile[grep("^[[:digit:]]", thefile[,2]),]


DiffCancer=unique(thefile[,1])
DiffCancer= DiffCancer[as.logical(nchar(DiffCancer))]

DiffCancer = DiffCancer[DiffCancer!="Total"]
DiffCancer = gsub(" [[:print:]]+$", "", DiffCancer)

Ncancer = length(DiffCancer)
Nages = length(unique(thefile[,2]))

Cancertype = thefile[,1] = rep(DiffCancer, rep(Nages, Ncancer))

# reformat the dataframe, with nice column names and       
mydata=data.frame(Cancertype = Cancertype, age=as.integer(substr(thefile[,2],1,2)), 
	Incidence=as.integer(gsub(",","", thefile[,4])), Rate= as.numeric(thefile[,5]),
	Population = as.numeric(thefile[,6]) )


theages = unique(mydata$age)
   
lambda = matrix(mydata$Incidence, nrow=length(theages), ncol=Ncancer, 
	dimnames = list(as.character(theages), DiffCancer)) 

CancerTotal = apply(lambda, 1, sum) 

lambda = lambda / matrix(CancerTotal, nrow=length(theages), ncol=Ncancer)

return(lambda)
}

lambdaDiffCVD=function(myfile,agegroup)   {
require(foreign)
thefile = read.table(file=myfile, sep=",", header=T, as.is=T, quote="\"")

# get rid of rows with totals in them

DiffCancer=c(unique(thefile[,1]),"others")
Ncancer = length(DiffCancer)

Nages = length(unique(thefile[,2]))

# reformat the dataframe, with nice column names and       
mydata=data.frame(Cancertype = thefile[,1] , age=as.integer(substr(thefile[,2],1,2)), 
	Incidence=as.integer(gsub(",","", thefile[,3])), Rate= as.numeric(thefile[,4]),
	Population = as.numeric(thefile[,5]) )


theages = unique(mydata$age)

others =as.numeric( mydata[mydata$Cancertype=="Total",]$Incidence -  mydata[mydata$Cancertype=="hf",]$Incidence-
mydata[mydata$Cancertype=="ihd",]$Incidence-
mydata[mydata$Cancertype=="stroke",]$Incidence )

otherspop = mydata[mydata$Cancertype=="Total",]$Population
othersrate = others/otherspop*100000
others = cbind(Cancertype=rep("others", length(theages)), age=as.numeric(theages), Incidence=others, Rate=othersrate, Population=otherspop)

mydata = rbind(mydata, others) 
lambda = matrix(as.numeric(mydata$Incidence), nrow=length(theages), ncol=Ncancer, 
	dimnames = list(as.character(theages), DiffCancer)) 

CancerTotal = lambda[,1]

lambda = lambda / matrix(CancerTotal, nrow=length(theages), ncol=Ncancer)

return(lambda)
}
